%% feature = extractDCTFeatures(coefarr, p, T)
% coefarr: 
%   M*N DCT-coefficient matrix. M = 8k, N = 8l (k,l are DCT-coef block
% numbers in each direction), the values in the coefarr must be previously
% thresholded.
% p: 4 1*N matrix (p.x1, p.x2, p.y1, p.y2 values between 1..8)
%   (p.x1, p.y1) -> (p.x2, p.y2)
% T:
%   threshold value
function feature = extractDCTFeatures(coefarr, p, T)
    %% Initiaize
    RNG = 2*T+1; % -T .. T
    len = length(p.x1);
    if (length(p.x2) ~= len) || (length(p.y1) ~= len) || (length(p.y2) ~= len)
        ME = MException('extractDCTFeatures:p', 'Matrix size does not match.');
        throw(ME);
    end
    value = zeros(RNG, RNG); % [-T .. T, -T .. T]
    [N, M] = size(coffarr);
    r1
    %% Extract feature block by block
    for k = 1:len
        for bi = 0:8:N-1
            for bj = 0:8:M-1
                v1 = coefarr(bi+p1i,bj+p1j);
                v2 = coefarr(bi+p2i,bj+p2j);
                value(v1+T,v2+T) = value(v1+T,v2+T)+1;
            end
        end
    end
    feature.value = value;
    feature.xrange = [-T, T];
    feature.yrange = [-T, T];
end
